function [X_mean, X_std] = Calc_quarterly_ave_NHz(Phi_NHz,ss_dist_NHz,X_NHz)
% computes the following:
% 1. E[(X(t)+X(t+1)+X(t+2))/3]=E[CondMean(t)] where
%    CondMean(t)=E[X(N(t),z(t))+X(N(t+1),z(t+1))+X(N(t+2),z(t+2))|N,z]
% 2. Var[(X(t)+X(t+1)+X(t+2))/3]=Var[Qave(t)] + E[CondVar(t)] where
%    CondVar(t)=Var[(X(N(t),z(t))+X(N(t+1),z(t+1))+X(N(t+2),z(t+2)))/3|N,z]

    CondMean = (X_NHz(:) + Phi_NHz*(X_NHz(:) + Phi_NHz*X_NHz(:)))/3;
    Cond2ndMoment_term1 = (X_NHz(:).^2 + Phi_NHz*(X_NHz(:).^2 + Phi_NHz*(X_NHz(:).^2)))/9;
    Cond2ndMoment_term2 = 2*(X_NHz(:).*(Phi_NHz*(X_NHz(:)+Phi_NHz*X_NHz(:))) ...
        + Phi_NHz*(X_NHz(:).*(Phi_NHz*X_NHz(:))))/9;
    CondVar = Cond2ndMoment_term1 + Cond2ndMoment_term2 - CondMean.^2;
    
    X_mean = sum(ss_dist_NHz(:).*CondMean);
    
    X_std = sqrt(max(sum(ss_dist_NHz(:).*CondVar) ...
        + sum(ss_dist_NHz(:).*(CondMean.^2)) - X_mean^2,0)); 
    % note: numerical rounding errors may result in negative variance when
    % variance is very close to 0.

end